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Abstract 

We study numerically classical 1-dimensional Hamiltonian lattices involving inter-particle long range 
interactions that decay with distance like l/r“, for a > 0. We demonstrate that although such systems 
are generally characterized by strong chaos, they exhibit an unexpectedly organized behavior when the 
exponent a < 1. This is shown by computing dynamical quantities such as the maximal Lyapunov 
exponent, which decreases as the number of degrees of freedom increases. We also discuss our numerical 
methods of symplectic integration implemented for the solution of the equations of motion together with 
their associated variational equations. The validity of our numerical simulations is estimated by showing 
that the total energy of the system is conserved within an accuracy of 4 digits (with integration step 
T = 0.02), even for as many as N = 8000 particles and integration times as long as 10® units. 


1 Introduction 

Hamiltonian systems describing 1-Dimensional (ID) particle chains characterized by interactions of various 
ranges constitute an active area of research with increasing interest due to their applicability in many 
scientihc fields. In particular, the relevance of long (vs. short) range interactions has been extensively 
studied and intensely debated in a wide variety of problems of statistical mechanics, mean field theories, 
active matter, dynamical networks, etc. regarding the various degrees of chaos involved in their time 
evolution. In statistical physics for example, the classical Boltzmann framework for the appropriate entropy 
functional at thermal equilibrium is not adequate for describing systems with long range interactions, as 
remarked already by J. W. Gibbs m- 

In the past 25 years, a great number of researchers [221 [H [m EH [241 [3] have shown that there exist 
long lasting quasi-stationary states (QSS) in a variety of physical and biological systems characterized as 
non-additive, i.e. that cannot be decomposed in entirely independent parts. In such cases, a different 
entropy functional (the so-called Tsallis entropy) appears to be more suitable for their thermodynamic 
description, while the associated probability distribution functions are of the g-Gaussian type with g > 1. 
This divergence from the classical Boltzmann-Gibbs Maxwellian distributions of g = 1, raises new questions 
regarding the statistical and dynamical behavior of such systems in the thermodynamic limit of very large 
N and total energy E, with E/N constant. 
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As we have demonstrated in recent publications mis], systems with long range interactions (LRI) have 
significant advantages, since they often exhibit a weaker ehaos than those with interactions only between 
nearest neighbors. For example, in active matter systems consisting of self-propelled particles (like birds) 
it has been observed that nonlocal communication acts as a counterbalance against external threats or 
attacks. This is the case, for instance with flocks of starlings described by the so-called topological model 
introduced by the STARFLAG group [3]. Their observations on groups of starlings in Rome revealed 
that synchronized movements are based on a fixed number of interacting neighbors, independent of the 
distance between them. Moreover, it was shown in [8] that the number n of interacting particles in the 
topological model is crucial for the coherence of the group: n needs to be large enough to overcome random 
perturbations as well as maintain cohesion. 

On the other hand, the study of LRI in the classical framework of dynamical systems is of great interest. 
Spatially localized oscillations called breathers, well-known as simple periodic solutions of lattices with 
nearest neighbor interactions, make their appearance also in Hamiltonian systems with LRI (see [13] for 
more details). Another interesting type of collective behavior occurs in the form of long-living QSS of 
the type mentioned above, while in a system of N coupled planar rotators that interact via long-range 
forces critical regions were found where these states appear m- Moreover, in a systematic study of 

the largest Lyapunov exponent showed that it decays as a power-law with N, making the system quasi- 
integrable in the thermodynamic limit. Such QSS were also studied in a generalized mean field system, 
where interactions decay with distance according to 1/r" jT] |T0| . 

More recently, the Fermi-Pasta-Ulam-/3 model (FPU-/3) with LRI was studied in (9] [I7| and opened a 
new branch of research in this field. In m for example, the authors explore instability regimes for the 
low frequency modes of an FPU-/3 chain in relation to the long-standing question of the relaxation times 
required to reach energy equipartition among all modes. 

In the present paper we focus on two topics: (i) First we extend recent studies of nearest-neighbor 
Hamiltonian lattices to analogous models involving LRI, and (ii) we employ numerical integration schemes 
to compute the tangent dynamics needed for the calculation of the largest Lyapunov exponent. The 
structure of our paper is as follows: In Section [2] we discuss in detail the transition to LRI for different 
boundary conditions, while in Section [3| we solve the associated variational equations for specific chaotic 
states. Subsection o gives the basic properties and dehnitions of symplectic integrators and discusses 
in Subsection 13.21 the so-called tangent map method. Finally, in Section HI we apply these techniques to 
two different Hamiltonian systems: the mean field model of planar rotors and the FPU-/3 chain, both with 
interactions modulated by the factor 1 /r“. We end with our conclusions in Section [5j 

2 Hamiltonian particle chains with long range interactions 

Let us consider a ID Hamiltonian chain with quadratic kinetic energy in the generalized momenta and 
a potential that depends purely on the generalized positions and contains nearest neighbor interactions 
together with an on-site potential. This system is described by the Hamiltonian 

■H = + + , (1) 

i i i 

where pi and Xi are canonical conjugate pairs of positions and momenta and boundary conditions are 
chosen arbitrarily. 

To convert the ID lattice ([T|) to an LRI system, all nearest neighbor interacting term^ in the potential, 
i.e the position differences Xj+i — should be replaced by the difference combinations Xi — Xj, i,j = 

^Or some of them. In [9] we considered LRI on the quartic potential of FPU-/3 model, while the quadratic one was left 
unchanged. 


2 




Figure 1: Lattices with nearest neighbor interactions (a —>• oo) on the left are converted into long range 
systems on the right. In our paper the strength of this interaction decays with distance as l/r“, for 
finite a > 0. Upper panels correspond to fixed or open boundary conditions and lower panels to periodic 
boundary conditions. 


1,..., A^, so that particles can form a fully connected graph (see Fig. [T]). Then, of course, one can apply 
interactions which depend on the topological distance and decay with distance as l/r“, with a > 0, as in 
[T]. Next, (dl) is converted to a Hamiltonian with LRI: 




+ E»'fe) + ^E 


1 M-^V{xi-Xj) 


( 2 ) 




where rjj defines the topological distance, i.e. the minimum connection between the particles i and 
j. N is the rescaling factor necessary for making the Hamiltonian ([2]) extensive, which means that the 
potential function V is proportional to the system siz^. In particular, when the chain has fixed or open 
boundary conditions, we set 


=\ i — j I and N 


1 

N 


N iV+1 

EE 


2=0 
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and when it has periodic boundary conditions we write: 


Xij = min{| i — j |,N— | i — j |} and N 


N/2-1 

2 E 


i=l 



(3) 


(4) 


The role of the parameter a which controls the range of interaction is essential. The two extreme cases 
are: (i) a —oo, where only interactions between first neighbors apply and (ii) a = 0, which corresponds 
to the mean field model (HMF), where all particles interact equally and independently of their distance, 
as in a fully connected graph. 

It appears that the implementation of LRI in a Hamiltonian particle chain strongly affects the system’s 
dynamics and statistics. The excitation of a single site results in an almost immediate reaction of all 

■^Otherwise the sum with long range terms would be approximately proportional to N and would converge faster to 
infinity than the other two sums of m in the thermodynamic limit. 
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the particles. The reaction time is significantly longer in simple chains. There is a number of physically 
interesting models whose dynamics can be extended to long range and thoroughly studied to provide 
representative examples for the LRI effect: (i) the Klein-Gordon chain for W{x) = ± and V{x) = 

in ([2|), (ii) a system consisting of coupled planar rotators with V{x) = 1 — cosx (and W{x) = 0) and 
(iii) the Fermi-Pasta-Ulam-/3 (FPU-/3) model with V{x) = + ^x^ in ([2]). Systems (ii) and (iii) are 

the ones treated in the present paper, while extensions of the KG model (i) will be dealt with in a future 
publication. 


3 Numerical integration of LRI systems and their tangent dynamics 

The equations of motion for the Hamiltonian system ([2]) with LRI read: 


Pi 


Pi 

dWjxi) 

dxi 





1 dV{xi — Xj) 
rfj dxi 


(5) 


The tangent dynamics (or variational equations) of a dynamical system is of great importance in 
uncovering the local properties of its solutions that lead to the estimation of chaotic indicators such as 
Lyapunov and GALI exponents [20]. In the case of a mechanical system the variational equations are 
defined on the tangent space of the configuration space, constituting a 4A^ dimensional phase space. In 
particular, the variational equations of the system ([2|) can be written in the following form: 

5y. \ ^ ( On In \ ( \ 

) V On J \ Sp J 


where each element of the matrix An is given by the second partial derivatives 



d^Wjxj) _ 1 d'^V{xi-Xk) 

1 1 d^V{xi-Xj) 

N dxidxj ’ 


if i = 3 

if i^j 


(7) 


calculated along the orbit (x(t),p(t)), while On, In are the zero and identity N x N matrices respectively. 

Now the equations of motion Q together with their variations (l6|) constitute a time-dependent Hamil¬ 
tonian system defined on the tangent space of the system ([2|) and expressed by: 


n 




1 y^ a^IT(xO 

2 ^ dx? ^ 


+ —E 

2N ^ 


1 d‘^V{xi — Xj) 
rf- dxidxj 


5xi5xj 


( 8 ) 


3.1 Symplectic integrators 

Symplectic integration schemes are designed to conserve the symplectic structure of a Hamiltonian system 
exactly, but not its total energy. During the computations, neither the orbit nor the Hamiltonian H 
itself are precisely followed. Even though these are not desirable properties for an integration method, 
symplectic integration still has great advantages compared with other methods [25]: (i) There exists 
an analytic autonomous Hamiltonian system H which is followed exactly |6| and provides the so-called 
backward error analysis, (ii) Energy errors are bounded and hence do not increase for exponentially long 
times. 
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For reasons of self-consistency and in order to help the reader who is not familiar with symplectic 
integration theory, we give here some notations in the spirit of Lie algebras [181 [25] . A symplectic integrator 
approximates the solution x{t) = of the system by splitting the kinetic and potential terms of the 

Hamiltonian, i.e. ^ where A := Dt is the kinetic and B := Dy the potential operator, 

which do not commute in general. If we denote by r the time step of the algorithm, the truncation in the 
Taylor expansion of determines the order of precision in terms of r. 

The standard 2nd order symmetric exponential splitting [2S], is: 

S2{t{A + B)) := + 0{t^) , 

with truncation errors of order r^. Furthermore, Yoshida has shown [25] that the symplectic integrators 
for n > 2 and n even are constructed by the composition of 2nd order ones. More specifically, he proved 
that the constants satisfying a 4th order symplectic scheme: 

Si{T{A + B)) = S2{KiT{A + B))-S2{n2T{A + B))-S2{KiT{A + B)) + 0{T^) , (9) 

are the solutions of the algebraic equations: 2ki + K 2 = 1, 2k\ + = 0. 


3.2 The tangent map method 


As already mentioned in SectionjHJ the variational equations of ([JD constitute a time-dependent Hamiltonian 
system and therefore admit a symplectic integration scheme (see m for more details). The extended 
Hamiltonian ([8]) is defined on the tangent space of the phase space solutions and consequently the operator 
Lh acts on the extended position momentum vector u = {x,p, 6x, 6p). In particular, Lh decomposes again 
into the kinetic-potential operators: 

Pk = Pk 

^Xk=Pk-\ + ^Xk 
6pk = dpk 

and 




Pk = -§f,-X+Pk 


— ^^k 


^Pk = - (Ef=i ■'^ + ^Pk 


through which we can now proceed to solve the variational equations of the problem. 


4 Applications to two important models 

4.1 Planar rotators and the mean field model 


The ID XY model is a chain of N planar rotators, describing a system of spins coupled by nearest neighbor 
interactions. The implementation of LRI to this model has already been studied in various interesting 
papers [Ml u uniisi n] in connection with its statistical and dynamical properties. The Hamiltonian of 
this system reads: 


n 


Y.p'i 


+ 


— E 

2N ^ 


1 — cos{'&i — 'dj) 




( 10 ) 
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where for a —)■ cx) the classical nearest-neighbor rotators are recovered, while for a = 0 (1101) reduces to the 
so-called HMF model. Regarding this model, Latora et al. m studied the a = 0 case, while Anteneodo 
and Tsallis [T] explored the more general a > 0 system. 

The variational equations for this model are easily derived from (1101) : 


where 


Pi 

Mi 

SPi 


Pi, 


1 sin(??j — -dj 

Wi 


^Pi 

3 


_L 

N 

1 cos(i?i-i?j) 
M r9. ’ 


if i = j 
if i^j 


( 11 ) 


The system (jlip is used for the calculation of the maximal Lyapunov exponent A. After a careful 
study of the behavior of the Lyapunov exponents in the HMF model, Latora et al. m established how 
they depend on the specihc energy e = U{N)/N and the system size. More specihcally, in a graph of A 
versus the specific energy, all data show a peak around e = 0.67, which weakly depends on the system size. 
Furthermore, A exhibits a decreasing level of chaos with increasing N as it decays as A r\j N In the 
same spirit, Anteneodo et al. in [T] also demonstrated a power law decay of A vs. N for the system (jlOp 
and any a < 1. 

As we discuss below, the FPU model under LRI exhibits a power-law decay of the maximal Lyapunov 
exponent with increasing N similar to the coupled rotators system. However, remarkably enough in the 
FPU case A does not decay as the specihc energy increases, but continues to grow as a power law. 


4.2 The FPU-/3 ID particle chain with LRI 

In the case of the FPU-/3 chain, LRI can be implemented in various ways since the system has a potential 
with two distinct nearest neighbor parts, a quadratic and a quartic one [9] 

N N N 

n=l n=l n=l 


Clearly, the most general approach would be to add LRI with different a exponents to each part of the 
potential as follows: 


n 


2 + 4JV 2-. + gjY Z-. ^“2 


U{N) 


(13) 


providing thus the system with two independent ranges, an ai controlling the quadratic part of the potential 
and an 0:2 controlling the quartic part. One can then focus on the computation of the maximal Lyapunov 
exponent A, as N increases at hxed specihc energy e = U{N)/N. To this end we need to integrate the 
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variational equations: 


with 


Xi 

= Pi, 


Pi 

1 Xi — Xj 

P (iCj - X 
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^Pi 
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/J-i _ 3j5 (xj— 

Ni N2 r“2 • 
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N2 


if i = j 
if i^j 


( 14 ) 


using the tangent map method of Section 13.21 with m = . The specific symplectic scheme we have 

employed is the Yoshida 4th order exponential splitting ([9]), while the method for calculating the Lyapunov 
exponent is based on Benettin et al. [5]. In Figl2]we plot the results of the model with LRI imposed only 
on the quartic potential, i.e. oi —>■ oo and a := a 2 > 0, so that only the long range interactions involved in 
(jl4p are nonlinear. Furthermore, in all computations of this section we apply fixed boundary conditions, 
while periodic boundaries give very similar results. The initial conditions correspond to all positions equal 
to zero and momenta drawn randomly from a uniform distribution. 

The evolution of the maximal Lyapunov exponent A is displayed in Fig. ETa) versus the system size 
N in double logarithmic scale for several nonlinear ranges a and energies that vary proportionally to N, 
keeping the specific energy s = 9 fixed. As pointed out in previous studies, we also find here that for a < 1 
the Lyapunov exponent decays as N increases by a power-law of the form X{N) ~ N~'^, k > 0 implying 
thus that the system tends to become less chaotic as the degrees of freedom grow. Instead for a > 1 
the Lyapunov exponents tend to stabilize at a constant value. The higher is the a value, the larger the 
A. Finally, for a = 10 we practically recover the classical FPU-,d model and as can be seen in Figl2](a) 
the symbols of these two cases superpose perfectly. Plotting then in Fig. [2](b) A as a function of time for 
N = 512,1024,2048,4096,8192, a 2 = 0.4 and ,5 = 10 we compute these decreasing values of A at t = 10® 
where they appear to stabilize for times greater than 10^. 

The challenging question, therefore, is whether A continues to decrease in the thermodynamic limit, 
where both N and U{N) increase indefinitely. This would be quite surprising since it would imply that 
the system would show signs of integrability in that limit! A more plausible alternative, of course, is that 
the dynamics tends to an “edge of chaos” state where it continues to be chaotic, albeit very weakly so, 
with a maximal Lyapunov exponent equal to zero. Since our computational capacity does not allow us to 
shed more light on this question, we shall have to postpone its investigation to a future analysis. 

We have also studied the dependence of the maximal Lyapunov exponent on the specific energy of the 
system. Recall that in subsection 14. 1 1 we mentioned that A has been found in the planar rotators model to 
peak at a specific e-value. To probing further this possibility, we studied in Figj3]the behavior of A versus 
e of the FPU-,5 model with a = 0.4. In contrast to the case of planar rotators, we observe that in the 
FPU system there is a monotonic increase of the maximal Lyapunov exponent with specific energy. In 
particular, A increases with the specific energy as A = 0.205e®'^^®. The same behavior is conjectured for 
all a values, since similar results have emerged for the classical FPU-,5 model as well. Nevertheless, the 
precise scaling laws depending on a will be studied in a future publication. 

Finally, as a test for the accuracy of our calculations, we compute the relative errors 


RE = log 


m) - m) 

m) 
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Figure 2: (a) The dependence of the maximal Lyapunov exponent on N, calculated at t = 10® for various 

ranges of nonlinear interaction a (e = 9, /3 = 10 and FPU-/3 means a = oo). (b) Evolution of the maximal 
Lyapunov exponent for a = 0.4, e = 9 and /3 = 10 for various system sizes N. 



Figure 3: Numerical calculation of the maximal Lyapunov exponent as a function of the specific energy 
for the system with N = 2048, a = 0.4 and /3 = 10. The red line corresponds to A = 0.205e®'^®®. 
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Figure 4: The relative errors of: (a) the Hamiltonian FPU-/3 with LRI and (b) the classical FPU-/3 
Hamiltonian. Different system sizes have been considered and same time step r = 0.02. 


where E(t) represents the total energy of the solutions as a function of time. Fig. Ill|a) shows the relative 
errors of the Hamiltonian (jl3l) for the LRI data of Figl2]^b), while FiglD^b) refers to relative errors in the 
FPU-/3 model with only nearest neighbor interactions. Our results verify in each case that the integrated 
Hamiltonian H remains within the same accuracy (nearly 4 significant digits) close to the desired Hamil¬ 
tonian H for a time step r = 0.02. In addition, the energy fluctuations become smaller and smaller in 
both cases when N increases. This is interesting because it implies that RE fluctuations are independent 
of the behavior of the maximal Lyapunov exponent, which is known to increase as a function of N for 
the classical nearest-neighbor FPU system (FiglU^b)), while it decreases for the FPU with LRI (FigHJ^a)). 
Nevertheless, we have no explanation why the energy fluctations are more pronounced in systems with 
long range interactions compared to nearest neighbor systems like FPU. 

5 Conclusions 

The study of the range of interactions in dynamical systems describing large numbers of particles is of great 
interest in many branches of science. In statistical physics, for instance, where the particles are modeled 
by nonlinear mechanical oscillators it is very important to understand the effect of mutual interactions on 
the dynamics and statistics of 1-dimensional Hamiltonian lattices in the so-called thermodynamic limit, 
where the number of particles N and their total energy E tend to infinity with e = E/N constant. For 
example, in a system of N coupled planar rotators that interact via long-range forces critical regions have 
been found in connection with long-living quasi-stationary states. On the other hand, such states were 
also studied in a generalized mean field system, where mutual interactions decay with distance according 
to l/r“, and it was found that the largest Lyapunov exponent decays as a power-law with N. 

In the present paper we have used advanced numerical techniques to extend previous studies of nearest- 
neighbor Hamiltonian lattices to analogous models involving long range interactions. To achieve this, we 
have employed symplectic integrator schemes and have computed the tangent dynamics needed to calculate 
the largest Lyapunov exponent of two different Hamiltonian systems: the mean field model of planar rotors 
and the FPU-/1 chain, both involving interactions whose range is modulated by the factor l/r“. We have 
thus been able to identify in these systems a continuous transition from strongly to weakly chaotic dynamics 
in the thermodynamic limit, as the parameter ol varies from infinity to zero. 
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We thus conclude that a = 1 appears to constitute a critical crossover value where a qualitative change 
occurs in the dynamics and statistics of such systems. As a chaotic index we have focused on the largest 
Lyapunov exponent A calculated efficiently and accurately using Yoshida’s 4th order symplectic integrator 
scheme. For both models treated in this paper we found that the Lyapunov exponent gives the same 
results: (i) For a > 1, A tends to stabilize at a positive value as N increases, (ii) for a < 1 A decreases 
with system size as for some positive constant k(q;) that depends on a. On the other hand, unlike 

the mean field rotator model, the A of the FPU-/3 chain increases monotonically as the specific energy 
£ = E/N IS increased. 

It is certainly somewhat intriguing that when long range interactions are applied only to the harmonic 
part of the potential no crossover is observed at a = 1, as the dynamics remains strongly chaotic and the 
statistics is always of the Boltzmann-Gibbs type [9]. On the other hand, all the studies carried out so 
far regarding long range interactions have been performed on Hamiltonian systems involving only mutual 
interactions and no on site potential. It would, therefore, be very interesting to study these cases also in a 
future publication to complete our understanding of the effect of the range of interactions in 1-dimensional 
Hamiltonian lattices. 
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